Ab-initio calculation of phonon dispersion curves: accelerating q point convergence 
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We present a scheme for the improved description of the long-range interatomic force constants in 
a more accurate way than the procedure which is commonly used within plane-wave based density- 
functional perturbation-theory calculations. Our scheme is based on the inclusion of a q point grid 
which is denser in a restricted area around the center of the Brillouin Zone than in the remaining 
parts, even though the method is not limited to an area around Y. We have tested the validity of our 
procedure in the case of high-pressure phases of bulk silicon considering the bet and sh structure. 
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I. INTRODUCTION 

The calculation of vibrational properties of solids, sur- 
faces, and nanocrystals plays an important role in the 
structural characterization of matter. Besides the pos- 
sibility to compare the vibrational frequencies with ex- 
perimental data, phonon modes can be used also to de- 
termine the structural stability of a system This 
requires a reliable description of the vibrational proper- 
ties of a system. Nowadays, ab-initio methods based on 
the density-functional perturbation theory (DFPT)£ or 
frozen-phonon techniques^ are commonly used for this 
purpose. With the latter the phonon frequencies can be 
obtained by calculating the dynamical matrix for each 
q point along the high-symmetry directions of the Bril- 
louin zone (BZ), which can quite cumbersome. In the 
former case the dynamical matrices are computed just 
on a (finite) grid of q points in the irreducible wedge of 
the Brillouin zone (IBZ) with subsequent Fourier inter- 
polation. Because of the computational effort of these 
calculations, the grid is often taken as small as possible 
resulting sometimes in an inaccurate description of the 
long-range force constants yielding erroneous frequencies 
especially in the low- frequency range near the T point. 
Nevertheless, also other frequencies in the Brillouin Zone 
can be affected. 

Considering the investigation of the thermodynamic 
stability of the system, an error in the description of the 
low- frequency area close to T plays just a minor role for 
the calculation of the free energy. However, for the in- 
vestigation of the Gruneisen parameters or the thermal 
expansion especially at low temperatures it is necessary 
to describe particularly these frequencies correctly, since 
their reciprocal value enter the formula.— Furthermore, if 
someone is interested in the determination of soft phonon 
modes, which point at the instability of structure, an ex- 
act description of the corresponding frequencies is neces- 
sary. In this context one has to mention that a wrong 
description of the long-range force constants may yield 
imaginary frequencies for a stable structure, which are 
interpreted as a result of a structural instability. 

A possibility to overcome this problem has been given 



by Gonze and Lee& by correcting the long-range dipolc- 
dipole interaction contribution to the force constants us- 
ing a term which yields the correct non analytic behavior 
in the limit of small q similar to the non analytic correc- 
tions yielding the LO-TO splitting. With this procedure, 
the description of phonons close to T is significantly im- 
proved and the required number of q points is reduced. 
However, this scheme is restricted to semiconductors and 
can not be employed in case of discrepancies at points 
close to the boundary of the BZ, too, at least this possi- 
bility has not been proven yet. 

We have taken another approach to overcome the prob- 
lem of the correct description of the long-range force con- 
stants within DFPT, which is related to the q-point con- 
vergence by introducing a mini-Brillouin Zone (mini-BZ). 
Inside of this mini-BZ the dynamical matrices are com- 
puted for a denser mesh of q points. The mini-BZ can 
be chosen around the center of the BZ but also at its 
boundary. These additional contributions are taken into 
account in the determination of the force constants. The 
validity of this procedure has been verified in the case of 
the high-pressure phases of bulk silicon: we have chosen 
the body-centered tetragonal structure (bet) correspond- 
ing to the /3-tin and the simple hexagonal (sh) structure 
corresponding to the sh phase. In both cases structures 
beyond the range of structural stability of the phase have 
been selected. Our method yields an improved descrip- 
tion of the phonon dispersion curves using less q points 
than the standard refinement of the grid. In detail, imag- 
inary frequencies for the bet structure near T standard 
procedures have been traced back to erroneous force con- 
stants whereas the soft phonon mode for the sh structure 
have been verified. 

This article is organized as follows: In Sect. [H] we de- 
scribe the theoretical framework of our method. Next, a 
short review of the technical details of our investigation 
is given (Sect. Hi ll). In Sect. IV] we apply our procedure 
to the bet (Sect. UVA) and the sh (Sect. IIVB]) structure 
of bulk silicon, where the results are compared and dis- 
cussed afterwards in Sect. IIV CI Finally we summarize 
and draw a conclusion (Sect. |V| . 
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II. THEORY 

In general, the phonon frequencies w(q) at a given q 
point in the BZ are obtained by diagonalizing the dynam- 
ical matrix D QQ / («;«', q), i.e., by solving the equation 
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where k, k! label the sublattices (the basis atoms), a, a' 
are the cartesian coordinates, and 13jv a is the 3-/Va x 3Na 
unitary matrix for TVa atoms in the cell. We focus on 
the Fourier-interpolation scheme to obtain the frequen- 
cies along the high-symmetry directions of the BZ. To 
this end, DFPT calculations of the dynamical matrices 
are performed for q points on a finite, regular grid of q 
points. These dynamical matrices are connected with the 
force constants matrices <1> QQ < ( ,) by a discrete Fourier 
transform (FT): 
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where R(Z, k) is the coordinate of the nth atom in the 
Ith cell. Using the common procedure within plane-wave 
codes one calculates the force constants by a FT from 
the dynamical matrices for the q points from the discrete 
and finite grid and subsequently obtains the dynamical 
matrices for any q point by FT of the force constant ma- 
trices. However, the drawback of this procedure is that 
the range of the forces is connected with the q points: 
the choice a lixl 2 xl 3 grid of q points leads to the inclu- 
sion of interatomic force constants between atoms within 
lixl 2 xl 3 cells. In other words, assuming a finer q point 
grid one can extend the range of the forces included in 
the force constants. Thus, usually convergence tests have 
to be performed comparing the phonon dispersion curves 
for calculations based on various grids. However, the 
number of q points is l\ ■ l 2 ■ I3 in the BZ (which can 
be reduced by symmetry), and therefore most investiga- 
tions are restricted to a smaller set of q points with a less 
accurate description of the long-range force constants. 

The long-range force constants affect particularly the 
low-frequency phonons close to the T point which might 
be wrongly described within this procedure. This prob- 
lem can be overcome for semiconductors by using the 
method described in Ref. [![. Another possibility is to 
use more q points in the area close to T which we have 
chosen here. 

Instead of taking increasingly grids we have taken a 
denser grid just in a small area around T (mini-BZ), and 
we assume a less dense grid outside this mini-BZ. Since 
the FT is based on a regular grid, we interpolate the 
missing dynamical matrices outside the mini-BZ by a FT 
from the force constants based on the coarser grid. In 
the following detailed description we neglect the indices 
k and a for simplicity, and we investigate the case l\ = 



h = h only. The case h ^ h ^ h follows analogously. 

Assuming an nxnxn grid, the force constants $ nn " 
are obtained by FT from the corresponding dynamical 
matrices B nnn calculated within DFPT: 
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From these force constants $ nTm the dynamical matrices 
for any q point in the BZ can be calculated by a back 
FT, therefore also for the q points on a finer grid, e.g., a 
Ixlxl grid with I > n. Thus, one gets 
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and also from these dynamical matrices D" lt (q) one can 
get again force constants $[" t by 



(5) 



which are in this case identical to In a next step 

we have performed calculations of the dynamical matri- 
ces D'"(q Gmini-BZ) for q points of the Ixlxl grid in- 
side the mini-BZ using the DFPT procedure. Taking 
these dynamical matrices and the interpolated ones out- 
side the mini-BZ, the improved force constants <&j" t can 
be calculated by 
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From these force constants $|" t the dynamical matrices 
for any q, in particular along the high-symmetry direc- 
tions of the BZ, are calculated in the standard way. In 
the limit of very fine grids one should achieve the same 
results as within the procedure of Gonze and Lee. 8 How- 
ever, this method is not restricted to a mini-BZ around 
the r point since the mini-BZ can be chosen arbitrarily. 
This has the advantage that also ranges of the disper- 
sion curves far away from the T point can be improved. 
Besides, the shape of the phonon density of states for 
particular spectral features can be inspected in detail. 

In order to apply this scheme we have modified a 
postprocessing routine of the QUANTUM ESPRESSO 
package^ in order to write out not only the frequencies 
but also the complete dynamical matrices after the FT. 
We have checked the numerical stability of the method 
by comparing the phonon dispersion curves based on the 
force constants $ 888 from Eq. ^ and $ 888 from Eq. (©, 
and we have found differences in the frequencies of less 
than 0.25 cm" . 

In the following we apply the procedure for testing 
purpose to the two silicon structures mentioned above 
and describe how to choose the mini-BZ and the required 
Ixlxl grid. 



III. METHOD 

All calculations have been carried out with the 
QUANTUM ESPRESSO package, 9 It is based on a 
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plane-wave pseudopotential approach to the density- 
functional theory (DFT)i±&±i For silicon we have em- 
ployed a norm-conserving pseudopotential generated fol- 
lowing the scheme suggested by v. Barth and Ca n 12 i 13 
The exchange-correlation energy is described within 
the local-density approximation (LDA) , 14 ' 15 We have 
used a kinetic-energy cutoff of 40 Ry and a 20x20x20 
Monkhorst-Pack mesh 16 together with a Methfessel- 
Paxton smearing^ using a width of 0.03 Ry to describe 
the electronic (metallic) ground state of the systems. 
The phonon frequencies have been calculated using the 
DFPT scheme±&±&2°. as implemented in the QUANTUM 
ESPRESSO package followed by a discrete Fourier Trans- 
form as described in Sect. HT1 

Both, the sh and bet structures have been investigated 
using a common body-centered orthorhombic cell (bco, 
lattice constants a ^ b ^ c) with two atoms at (0, 0, 0) 
and at (0,0.56, Ac) in the unit cell. The symmetry of 
bet requires a = b and A = 0.25 whereas the symmetry 
of sh yields b — \/3c and A = 0.5. In fact, for sh we 
use a biatomic supercell although the structure of the 
sh phase can be described with just one atom in the sh 
unit cell. However, using the bco cell we have access to 
soft modes corresponding to the doubling of the unit cell. 
For details of the choice of the cell see Ref. In this 
work, we have relaxed the ground-state geometry of the 
structure for a volume fixed to 184 a| for both sh and 
bet. The equilibrium lattice constants are c/a = 0.5489 
for bet and c/a = 0.5338 and b/a = 0.9230 for sh. The 
error with respect to the ideal b/c ratio of sh is negligible. 



IV. RESULTS 

For the application of our method we have chosen the 
bet structure of silicon at V = 184 a% which is a vol- 
ume beyond the stability range of the corresponding (3- 
tin phase. For this structure we have found a phonon 
instability along the T-X direction of the bet B Z 21 i 22 
which is equivalent to the T-T direction of the bco BZ 
(see Ref. [4(). This phonon instability turned out not to 
be physical. The phonons in the mentioned articles had 
been calculated using a 4x4x4 Monkhorst-Pack mesh, 
which was slightly insufficient to describe the frequen- 
cies in this region of the BZ properly, since calculations 
within DFPT of dynamical matrices at q points near T 
have yielded only real frequencies. Because the phase 
space of the numerically soft modes was negligibly small, 
the imaginary frequencies did not affect the results for 
the free energy. However, we use this case to check the 
validity of the method described in Sect. [IT] by the use of 
the mini-BZ based on a 4x4x4 grid outside and a 8x8x8 
one inside the mini-BZ. A comparison of the phonon dis- 
persion curve using the 4x4x4 grid plus the mini-BZ 
with to the one based on 8x8x8 mesh gives an estimate 
of the errors using our scheme. The bet structure is also 
taken to exemplify the choice of the extent of the mini- 
BZ. Then we apply our scheme to go beyond the 8x8x8 



[100] [010] [001] [110] [101] [011] [111] 



400 
















A 
\l' 


■ 


If 
\ r 














300 






Jr 


JjTw 










AVI 

Pi 


200 














I'A'A 


i'rW 
tr A\ 


w 


100 


~ Vtffe 

v» 




ft 
1 


s * 

\ 

\ 




iff VM 


V \ 




a' 
at 







Mi 


--2x2x2 
4x4x4 


\1 


Iff/ wll 






f BCT 






— 8x8x8 







x ruxHTTrRrsrwx 



FIG. 1: (Color online) Phonon dispersion curves for the bet 
structure at V = 184 a|. Dispersion curves obtained by FT 
using a 2x2x2 mesh (dashed lines), 4x4x4 mesh (dotted 
lines), and a 8x8x8 mesh (solid lines). 

mesh for the final results. 

As a second example we have chosen the sh structure 
of silicon at V = 184 a? B which is also a volume beyond 
the stability range of the corresponding phase, here the 
sh phase. Choosing a biatomic supercell, a soft phonon 
mode has been found at the T and the S point of the 
bco BZ. 4 Both points are equivalent to the T point of 
the monatomic sh unit cell. The finding of a soft phonon 
mode at the T point is in accordance with other reported 
resultsi 23 i 24 ' 25 In fact, the softening at the BZ-boundary 
point S refer to a doubling of the sh unit cell. This is in 
agreement with the distortion of the biatomic supercell 
which contains two monatomic sh cells. Thus, this soft 
phonon mode is of physical origin. 

Both case studies, the one with an unphysical but nu- 
merically imaginary frequencies and the one with the 
physically correct phonon instability are described in the 
following, and they are finally compared and discussed. 

A. Application to the bet structure of silicon 

1. Standard procedure: Increasing mesh size 

For the bet structure at V = 184 a\ corresponding 
to a pressure of 133 kbar the frequencies along the high- 
symmetry directions of the bco BZ have been calculated 
using a 2x2x2, a 4x4x4, and a 8x8x8 grid. The results 
are shown in Fig. [TJ Since the bet structure has a higher 
symmetry than the bco structure, the T-X and the T-R 
directions (bco) are equivalent to the T-U-X and the T-S 
directions (bet), respectively. The equivalent directions 
are shown mainly for completeness. The T-T direction 
with T(i, i, 0) (coordinates in units of reciprocal lattice 
vectors; in the following we will assume the points in the 
BZ always in units of reciprocal lattice vectors without 
mentioning it explicitly) is of particular interest because 



4 



50 - 
40 - 



E 30 - 



< 20 - . 



[100] [010] [001] [110] [101] [011] [111] 

i i ~r 



• 12x2x2-4x4x41 
o 14x4x4-8x8x81 

♦ 18x8x8-16x16x161 



10 - i 



Ik 



• o oj 



BCT- 



x ruxHrTTRrsrwx 

FIG. 2: (Color online) Differences Aui between the calculated 
and Fourier interpolated points along the high-symmetry di- 
rections for various meshes: Auj = I1J222 — ^444 PT | are drawn 
with solid symbols, Auj — \uj444 — u>^ FT \ with open symbols, 
and Auj — |w 8 gg — Wieil'iel are with crosses (see text). 



a soft phonon mode appears close to T using the 4x4x4 
mesh as visible in Fig. [T] Note, that this softening does 
not appear for the 2x2x2 grid, which is obviously insuf- 
ficient to describe the phonon dispersion correctly. In- 
creasing the mesh size, the softening remains for a 8 x 8 x 8 
grid, however, to a minor extent. From the dispersion 
curves it is difficult to decide whether convergence with 
respect to q has been achieved for the 4x4x4 or not, since 
the the frequencies at the mesh points of the 8x8x8 grid 
are on top of the Fourier-interpolated 4x4x4 dispersion 
curves. Only for the low-frequency mode along the T-T 
direction some of the interpolated frequencies are imagi- 
nary but the 8x8x8 points yield just real values for the 
frequencies. Besides, the shape of the 8x8x8 Fourier- 
interpolated phonon curves shows just minor deviations 
from the 4x4x4 ones. Inspecting the frequencies at mesh 
points of the 16x16x16 grid, the calculated frequencies 
are nearly indistinguishable from the interpolated disper- 
sion curves except along the critical T-T direction where 
all calculated points yield real frequencies, while a part 
of the 8x8x8 Fourier-interpolated phonon dispersion are 
imaginary. 

However, such a detailed study is not always possi- 
ble for every system. In our case, the 2x2x2 mesh re- 
quired the calculation of dynamical matrices at 4 q points 
in the IBZ, the 4x4x4 mesh at 13 q points, and the 
8x8x8 mesh at 59 q points. Since the dispersion curves 
do not change significantly assuming a grid denser than 
the 4x4x4 except along the T-T direction and there only 
in the region close to T it is not necessary to calculate 
all the dynamical matrices on a 8x8x8 or a 16x16x16 
grid. Ultimately, the unphysical instability should be 
lifted. This can be achieved by applying our method de- 
scribed in Sect. [TT] using a mini-BZ around T. The extent 
of the mini-BZ can be determined by inspecting the dif- 
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FIG. 3: (Color online) Phonon dispersion curves along F-T 
using a 4x4x4 and a 8x8x8 mesh together with interpolated 
curves (4x4x4m) based on various mini BZs as denoted in the 
inset. Imaginary frequencies are drawn along the negative 
frequency axis. The cutoff for the q points is in units of 
reciprocal lattice vectors for qi with i = x,y,z (see text). 



ferences between the frequencies derived from the DFPT 
dynamical matrices w DFPT and the Fourier-interpolated 
ones w lnt . For this purpose we have drawn in Fig. [5] the 
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using an Ixlxl and an nxnxn grid with I > n for the 
frequencies presented in Fig. [TJ As mentioned above, the 
differences Alu for the 2x2x2 and the 4x4x4 mesh are 
quite large. All differences decrease significantly using 
a finer grid, except for the points near T along the T- 
T direction. Thus, the application of our method for a 
mini-BZ around T promises an improvement of the re- 
sults especially in the range of T-T. 



2. Approval of the present Mini-BZ method 

For a first test we want to apply our method using a 
4x4x4 grid outside the mini-BZ and an 8x8x8 inside. 
The scope of this test is to reproduce the 8x8x8 curves 
(inclusive the numerically instable mode) using a 4x4x4 
mesh together with a Mini-BZ, since a phonon dispersion 
curve using a full 8x8x8 grid exists as a reference. 

Inspecting Fig. ® the largest deviation of the 8x8x8 
mesh from the 4x4x4 mesh (open symbols in Fig. \2§ of 
Aw w 36 cm' 1 is found for the point (|, |, 0) along T-T, 
but also for the point (|,|,0) the differences between 
the meshes are in the order of Auj ss 15 cm -1 . Besides, 
along X-r-U-X the differences are also remarkable for q 
points with components qi < 4 (i = x,y,z). Thus, the 
selection of q points up to qi < | would be a promising 
choice for the mini-BZ. 

In the following we have used various mini-BZs up to 
Qi < \ (i = x ,y, z ) for the interpolation (denoted as 
4x4x4m). The results for the low- frequency range of 
the phonon dispersion curve along the T-T direction are 
shown in Fig. [3] in comparison with the results based on 
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FIG. 4: (color online) Comparison of the dispersion obtained 
from a 8x8x8 grid (dotted line), using the qi < 0.375 mini BZ 
in addition to the 4x4x4 grid (4x4x4m, dashed line), and 
using 16x16x16 points in the qi < 0.25 mini BZ based on 
force constants from the 4x4x4m grid (8x8x8m, solid line), 
see text. Imaginary frequencies are drawn along the negative 
frequency axis. 

4x4x4 and 8x8x8 grids. Note: the choice of qi < 1 
would make the 4x4x4m grid identical to the 8x8x8 
grid. The curves for the mini-BZ with qi < | and 
qi < i are both very close to the curve using the full 
8x8x8 mesh. Therefore, the choice of < | for the 
mini-BZ which has been already assumed from Fig. [21 is 
confirmed. Also the high-frequency region of the disper- 
sion using an 8x8x8 grid is reproduced very well with 
this 4x4x4m grid (see Fig. [4j. Only minor deviations 
appear at regions more distant from T resulting from 
the unresolved Aw = 6.93 cm -1 along X-H-r. However, 
the general improvement of the accuracy of the phonon- 
dispersion curve is remarkable. Note, that with this mini- 
BZ only 15 q points of the 8x8x8 have been necessary 
in addition to the 13 q points of the 4x4x4 mesh, which 
are much less than the 59 q points of the full 8x8x8 grid. 



3. Application of the Mini-BZ to denser grids 

After verifying the validity of our method we want to 
go beyond the 8x8x8 mesh because of the remaining 
unphysical softening along T-T, whereas the calculated 
frequencies at 16x16x16 mesh points along the high- 
symmetry directions show only real values. Inspecting 
the differences Aw in Fig. [2] again there is just a ma- 
jor difference of ~ 15cm -1 along T-T for q= (-^, ^,0), 
whereas the other differences are tiny. Since along the 
T-T direction at q= (^,^,0) there is a crucial differ- 
ence of 3.3 cm -1 , which might be important for resolving 
the numerical soft mode, we have tested Mini-BZs with q 
points up to qi < j yielding additional 29 q points. The 
results are denoted as 8x8x8m. With this mini-BZ the 
phonon frequencies are described accurately as shown in 
Fig. Q] and the softening along T-T has been lifted. The 
curves for 8x8x8m and 4x4x4m match nearly exactly 



indicating that convergence has been achieved. Note, 
that in this case the mini-BZ using 16x16x16 points is 
applied on top of the mini-BZ using 8x8x8 points in 
addition to the 4x4x4 mesh. In fact, we achieve conver- 
gence for 8 x 8 x 8m since the remaining differences Aw are 
less than 1.25 cm -1 . It has to be mentioned that further 
improvement could be achieved using a mini-BZ close to 
the points T, R, and S, which is also possible within our 
scheme. 

In summary, we have been able to obtain converged 
phonon dispersion curves using in total 57 q points in 
the IBZ which is around one sixth of the 349 q points 
which arc required for a complete 16x16x16 mesh. In 
this way, the convergence is accelerated significantly. 

B. Application to the sh structure of silicon 

Similarly to bet we have investigated the sh structure 
of bulk silicon at a volume of V = 184 a^ which accords 
here to a pressure of 107 kbar, again beyond the range of 
stability of the corresponding sh phase. 

1. Standard procedure: Increasing mesh size 

First, we compare the phonon-dispersion curves using 
the 2x2x2, the 4x4x4, and the 8x8x8 grid and the 
differences Aw (see Eq.©) in Fig. [5] Because of the 
lower symmetry of the sh structure, more q points for 
each considered grid had to be calculated within DFPT: 
5 points for 2x2x2, 18 for 4x4x4, 95 for 8x8x8 and 621 
for 16x16x16. Inspecting Fig. O there are remarkable 
differences visible in the results based on a 2x2x2 and 
a 4x4x4 mesh and between the ones based on a 4x4x4 
and a 8x8x8 mesh. This is reproduced in the graph of 
Aw. However, the variations along the T-T direction are 
rather small and the the imaginary frequencies nearly do 
not change the extension. In addition to the T point 
there are significant differences around the X and the S 
points. Since the T point and the S point are in this case 
equivalent, an improvement at T will yield an improve- 
ment at S. Now we focus on the area around the T point. 



2. Second test of the present method 

Looking at the differences Aw between the 4x4x4 and 
the 8x8x8 mesh, the choice of qi < | as for bet is not rea- 
sonable for sh since the the deviation of Auj < 18.57 cm -1 
along the X-H~r direction at q= (0,0, |) cannot be re- 
duced with a mini-BZ around T. Nevertheless, there 
are differences of Aw in the order of 15 cm -1 close to 
r at q with qi < j which can be resolved. However, 
there are also significant deviations of sa 12 cm" 1 close 
to the R and S points. Thus we have chosen q t < i. 
With this mini-BZ in addition to the 18 q points of the 
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FIG. 5: (color online) Phonon dispersion curves for the sh 
structure at V = 184 a|. Upper panel: curves obtained from a 
2x2x2 mesh (dashed lines) , a 4 x 4 x 4 mesh (dotted lines) , and 
a 8x8x8 mesh (solid lines). Imaginary frequencies are drawn 
along the negative frequency axis. Lower panel: differences 
Aui for sh as in Fig. [2] 
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FIG. 6: (color online) Comparison of the dispersion obtained 
from a 8x8x8 grid (dotted line), using the qi < 0.5 mini BZ 
for 8x8x8 points in addition to the 4x4x4 grid (4x4x4m, 
dashed line), and using 16x16x16 points in the qi < 0.125 
mini BZ in addition to 4x4x4m (see text). Imaginary fre- 
quencies are drawn along the negative frequency axis. 



the 4x4x4m and the 8x8x8m grid. 

Considering the numerical effort, we have used 18 dy- 
namical matrices of the 4x4x4 mesh, additional 40 of the 
8x8x8 one, and furthermore 10 of the 16x16x16 grid, 
in total 68, which is much less than the 621 dynamical 
matrices required for a complete 16x16x16 grid. Also 
in this case the numerical effort has been reduced dras- 
tically. 



4x4x4 in further 40 dynamical matrices are necessary. A 
comparison between the phonon dispersion curves based 
on the full 8x8x8 mesh and the one using the mini-BZ 
(4x4x4m) is presented in Fig. [5] The agreement is ac- 
ceptable, only the overbending close to the X point is 
not described correctly due to the choice of the mini-BZ 
around T. Using an additional mini-BZ close to X would 
solve this problem. However, the scope here is a check 
of the method for a mini-BZ around T analogously to 
the bet case (Sect. IIV Al) . and this region is described 
excellently. 



3. Application of the Mini-BZ to go beyond the 8x 8x 8 grid 

Next, we want to go beyond the 8x8x8 mesh, again 
focussing on the region around the center of the BZ. The 
differences Aw between the 8x8x8 and the 16x16x16 
mesh show significant deviations for q^ < ^ , with a max- 
imum of Alo < 6.13 cm -1 at the S point, but also for the 
T-X direction along [100]. In order to reduce these dis- 
crepancies and for describing especially the range of the 
phonon softening correctly, we have chosen qi < | for the 
mini-BZ by including additional 10 dynamical matrices 
for the calculation of the force constants. In this way, 
the differences 6.05 cm -1 are eliminated along T-T. The 
resulting phonon-dispersion curve is presented in Fig. [5] 
One can notice the improvement of the results based on 



C. Discussion 

For both systems, the bet and the sh structure of 
bulk silicon, we have been able to apply successfully 
our scheme and we have obtained converged phonon- 
dispersion curves using less q points than the correspond- 
ing full mesh. Comparing the results of the bet and the 
sh structure, one can see that the convergence of the sh 
structure with respect to the q points is slower than the 
one of the bet structure which results in a larger choice 
of the mini-BZ. In particular, for bet the critical area 
with large variations with respect to the choice of the q 
mesh is around the center of the BZ whereas for sh it is 
around X. Since the ground state of both structures had 
been calculated using the same convergence parameters 
and the same unit cell at the same volume, the different 
speed of the q convergence is not due to different conver- 
gence parameters for the ground state. Therefore, one 
can not estimate the size of the q point grid which is re- 
quired for convergence from one structure to another one 
even using the same cell. Thus, the q point convergence 
has to be tested for any structure separately. However, 
it was possible to confirm the soft phonon mode at T for 
the sh structure, whereas the one of the bet structure 
has disappeared including enough q points around the 
center of the BZ. Hence, the latter phonon instability is 
due to an inaccurate description of the long-range force 
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constants and thus it has only numerical origin. It would 
be possible to reduce the remaining discrepancies for sh 
around the X point using our method by applying an ad- 
ditional mini-BZ. However, this is beyond the scope of 
this article. 



V. CONCLUSIONS 

We have presented a scheme within standard DFPT 
calculations for the accurate calculation of phonon dis- 
persion curves by improving the interatomic force con- 
stants which can be applied for semiconducting systems 
as well as for metallic ones. Especially the long-range 
contribution to the force constants can be described suc- 
cessfully. This scheme is based on the inclusion of a 
denser q-point mesh in some part of the BZ (mini-BZ) 
and a wider one outside. The method has been applied 
successfully to the bet and the sh structure of bulk sili- 
con, where the origin of a phonon instabilities has been 
discussed. In detail, for the bet structure the soft phonon 
mode has been traced back to an inaccurate description 
of the (long-range) force constants and the imaginary fre- 



quencies become real applying our procedure till conver- 
gence. In the case of the sh structure the soft phonon 
mode has been confirmed. For both cases, the number 
of required dynamical matrices has been reduced drasti- 
cally. Whereas here the mini-BZ has been chosen around 
the center of the BZ, our scheme allows also a different 
choice. In this way, also features at the boundary of the 
BZ can be described more accurately. The use of our 
method can improve results not only for phonon disper- 
sion curves but also for related quantities like Griineisen 
parameters, thermal expansion, and the free energy. 
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